# This code replicates the results for 
# Study 1 of Inclusionary & Exclusionary Preferences 

# Version: September 19, 2025 
# Author: Marika Landau-Wells (mlw@berkeley.edu)

# This file contains code for replicating:
# Main Text: Figure 2 & Figure 3
# Supplement: Tables C6-C7, D1-D14, 
# Modeling of raw fMRI data was conducted in neuroimaging-specific software (e.g., fsl)
# and is not reported here.  
# Group-level brain maps are on neurovault: 
# https://identifiers.org/neurovault.collection:13022

# For ease of interpretation, hypothesis labels follow the 
# article's format instead of the pre-registration
# H1a = pre-reg H-E1, H1b = pre-reg H-O1, H1c = pre-reg H-G1, 
# H2a = not pre-registered, H2b = pre-reg H-G1,
# H3a = pre-reg H-O2, H3b = pre-reg H-G2
# Additional hypotheses reported in Appendix E follow pre-reg labels

# Groups are always identified as: 
# Opponents = Group 1 
# Endorsers = Group 2

# Originally created with R version 4.4.3 (runs as of April 9, 2025)
# Running groundhog will ensure the correct versions of all packages are loaded
# for replication purposes

library(groundhog)

pkgs = c("tidyverse", "stringr", "ggplot2", "DescTools", "WRS2")
groundhog.library(pkgs, "2025-04-09")

# To load independently, use these versions:
#library(tidyverse) #v 2.0.0
#library(stringr) #v 1.5.1
#library(ggplot2) #v 3.5.1
#library(DescTools) #v 0.99.60
#library(WRS2) #v 1.1-6

# Set the working directory
setwd()

# Contents
# 0. Helper Functions
# 1. Import & Tidy Data
# 2. Define Search Space Datasets 
# 3. Normality Checks 
# 4. M1 Testing (H1a, H1b, H1c)
# 5. M2 Testing (H2a, H2b)
# 6. M3 Testing (H3a, H3b)
# 7. Additional Pre-registered Hypothesis Tests 
# 8. Figures (Fig. 2a & Fig. 3a)
# 9. Tables (C6-C7, D1-D14, F1) 

# Clear the environment
gc()
rm(list = ls())


################################
# 0. Helper Functions ####
## Within-Group robust, one-sided bootstrap-t method (alternative = "g") ####
# This is an adaptation of the trimcibt (bootstrap-t method) function in WRS2
# The objective is to take in a vector of difference scores and get p-value
# for Diff > 0 using the trimmed mean
trdiff1sp <- function(x,
                      nv = 0,
                      tr = 0.2,
                      alpha = 0.05,
                      nboot = 1000) {
  set.seed(1234)
  x <- as.numeric(x)
  test <- (mean(x, tr) - nv) / WRS2::trimse(x, tr)
  data <- matrix(sample(x, size = length(x) * nboot, replace = TRUE), nrow = nboot) - mean(x, tr)
  top <- apply(data, 1, mean, tr)
  bot <- apply(data, 1, WRS2::trimse, tr)
  tval <- sort(top / bot)
  icrit <- round((1 - alpha) * nboot)
  ibot <- round(alpha * nboot / 2)
  itop <- nboot - ibot
  p.value <- (sum(test <= tval)) / nboot
  
  result <- list(
    estimate = mean(x, tr),
    test.stat = test,
    tr = tr,
    p.value = p.value,
    n = length(x),
    alpha = alpha
  )
  
  class(result) <- "trimcibt"
  result
}



## Between-Group standard t-test for M2 data ####
grpttest.tom <- function(x, z){
  #x is the dataframe; z is the alternative
  #note z is specified to aid in directional hypothesis testing
  g1 <- filter(x, Group == 1) 
  g2 <- filter(x, Group == 2)
  rois <- unique(x$ROI)
  out.list <- list()
  for (i in 1:length(rois)){
    g1vox <- filter(g1, ROI == rois[i]) 
    g2vox <- filter(g2, ROI == rois[i]) 
    out.list[[i]] <- t.test(g1vox$mean_extracted_betas, 
                            g2vox$mean_extracted_betas, 
                            alternative = z, paired = F)
  }
  return(out.list)
}


## Between-Group Yuen's Robust T-test for M3 data #### 
# Note: Test direction should be set within the function
grp.robust <- function(x, y){
  #x is the dataframe; y is the contrast
  g1 <- filter(x, Group == 1, contrast == y) 
  g2 <- filter(x, Group == 2, contrast == y)
  rois <- unique(x$roi)
  out.list <- list()
  for (i in 1:length(rois)){
    g1vox <- filter(g1, roi == rois[i]) 
    g2vox <- filter(g2, roi == rois[i]) 
    out.list[[i]] <- DescTools::YuenTTest(g1vox$mean_topvoxels, 
                                          g2vox$mean_topvoxels, 
                                          alternative = "greater")
  }
  return(out.list)
}



## Between-Group standard t-test for M3 data ####
grpttest <- function(x, y, z){
  #x is the dataframe; y is the contrast; z is the alternative
  #note z is specified to aid in directional hypothesis testing
  g1 <- filter(x, Group == 1, contrast == y) 
  g2 <- filter(x, Group == 2, contrast == y)
  rois <- unique(x$roi)
  out.list <- list()
  for (i in 1:length(rois)){
    g1vox <- filter(g1, roi == rois[i]) 
    g2vox <- filter(g2, roi == rois[i]) 
    out.list[[i]] <- t.test(g1vox$mean_topvoxels, 
                            g2vox$mean_topvoxels, 
                            alternative = z, paired = F)
  }
  return(out.list)
}



################################

################################
# 1. Import & Tidy data ####

## Subject Data ####
s <- read.csv("Subject_Data_Shared.csv", 
              stringsAsFactors = T)
# Note: SUB_ID doesn't match neuro data, so need to create new variable to join dfs
slabel <- str_remove_all(s$SUB_ID,"_")
s$participantID <- paste("sub-",slabel,sep="")
#str(s)

# Reduce to matched subjects and columns for joining 
s.use <- filter(s, match == 1) %>%
  select(participantID, Group) 

## Multivariate data (for M1) ####
m1 <- read.csv("M1_Multivariate.csv", stringsAsFactors = T)

# Join with sub data
m1 <- inner_join(x=m1, y=s.use, by = "participantID")
#length(unique(d.m1$participantID)) #n = 42

# Pivot wide to create averages
m1.wide <- pivot_wider(m1,
                       id_cols = c(participantID, Group, roi),
                       names_from = c(contrast1, contrast2),
                       values_from = avg_z)
# Rename for ease
m1.wide <- rename(m1.wide,
                  NN = con_6_cgn_con_6_cgn,
                  ND = con_6_cgn_con_7_cgd,
                  NF = con_6_cgn_con_8_cgf,
                  NT = con_6_cgn_con_9_tgn,
                  DN = con_7_cgd_con_6_cgn,
                  DD = con_7_cgd_con_7_cgd,
                  DF = con_7_cgd_con_8_cgf,
                  DT = con_7_cgd_con_9_tgn,
                  FN = con_8_cgf_con_6_cgn,
                  FD = con_8_cgf_con_7_cgd,
                  FF = con_8_cgf_con_8_cgf,
                  FT = con_8_cgf_con_9_tgn,
                  TN = con_9_tgn_con_6_cgn,
                  TD = con_9_tgn_con_7_cgd,
                  TF = con_9_tgn_con_8_cgf,
                  TT = con_9_tgn_con_9_tgn)
# Create avg values across paired correlations
m1.wide <- m1.wide %>% 
  mutate(TN_avg = (NT + TN)/2,
         DN_avg = (ND + DN)/2,
         FN_avg = (NF + FN)/2,
         TD_avg = (DT + TD)/2,
         TF_avg = (FT + TF)/2,
         DF_avg = (FD + DF)/2,
         DDFF_avg = (DD + FF)/2,
         DDNN_avg = (DD + NN)/2,
         DDTT_avg = (DD + TT)/2,
         FFNN_avg = (FF + NN)/2,
         FFTT_avg = (FF + TT)/2,
         TTNN_avg = (TT + NN)/2)
glimpse(m1.wide)

# Haxby et al. (2001) difference scores:
m1.wide$wb_DF <- m1.wide$DDFF_avg - m1.wide$DF_avg
m1.wide$wb_DN <- m1.wide$DDNN_avg - m1.wide$DN_avg
m1.wide$wb_FN <- m1.wide$FFNN_avg - m1.wide$FN_avg
m1.wide$wb_TN <- m1.wide$TTNN_avg - m1.wide$TN_avg
m1.wide$wb_TD <- m1.wide$DDTT_avg - m1.wide$TD_avg
m1.wide$wb_TF <- m1.wide$FFTT_avg - m1.wide$TF_avg

# Back Transform to the Values used in Hypothesis tests
m1.wide$TD_r <- DescTools::FisherZInv(m1.wide$TD_avg)
m1.wide$TF_r <- DescTools::FisherZInv(m1.wide$TF_avg)
m1.wide$TN_r <- DescTools::FisherZInv(m1.wide$TN_avg)
m1.wide$DN_r <- DescTools::FisherZInv(m1.wide$DN_avg)
m1.wide$FN_r <- DescTools::FisherZInv(m1.wide$FN_avg)
m1.wide$DF_r <- DescTools::FisherZInv(m1.wide$DF_avg)



## Univariate data, Theory of Mind Network (for M2) ####
m2 <- read.csv("M2_Univariate.csv", stringsAsFactors = T)

# Subs (1 invalid Theory of Mind Localizer)
tom.subs <- unique(m2$participantID)
length(tom.subs) #41

# Recode Group to match other datasets
m2$Group <- ifelse(m2$group == "OPPOSE", 1,2)
m2 <- select(m2, -group)

# Add labels to ROI factor
m2$ROI <- factor(m2$ROI,
                 levels = c("DMPFC_binary", "LTPJ_binary",
                            "MMPFC_binary", "PC_binary",
                            "RSTS_binary", "RTPJ_binary",
                            "VMPFC_binary"),
                 labels = c("dmPFC", "l. TPJ",
                            "mmPFC", "Precuneus",
                            "r. STS","r. TPJ","vmPFC"))

## Univariate data, CBMA clusters (for M3 + other hypotheses) ####
m3 <- read.csv("M3_Univariate.csv", stringsAsFactors = T)

# Generate 1 value per subject per roi per condition by averaging folds
m3.use <- select(m3, 
                 participantID, roi, contrast, mean_topvoxels_leftout_cope)
m3.use <- m3.use %>% 
  group_by(participantID, roi, contrast) %>%
  summarise(mean_topvoxels = mean(mean_topvoxels_leftout_cope))

#42*23*7 #6762 = correct number of observations (all matched subs x rois x contrasts)

# Tidy the factors
m3.use$contrast <- factor(m3.use$contrast, 
                          levels=c("con_1_tgn-gt-cgn", "con_2_cgd-gt-cgn", "con_3_cgf-gt-cgn",
                                   "con_4_cgf-gt-cgd", "con_5_cgd-gt-cgf",  
                                   "con_10_tgn-gt-cgd", "con_11_tgn-gt-cgf"),
                           labels=c("tgn-gt-cgn", "cgd-gt-cgn", "cgf-gt-cgn",
                                    "cgf-gt-cgd", "cgd-gt-cgf", "tgn-gt-cgd", "tgn-gt-cgf"))
m3.use$roi <- factor(m3.use$roi, 
                     levels=c("emo_reg_cluster_lCingG", 
                              "emo_reg_cluster_lIFG", "emo_reg_cluster_lMFG",
                              "emo_reg_cluster_lMTG", "emo_reg_cluster_lSFG", 
                              "emo_reg_cluster_lTG", "emo_reg_cluster_rInsIFG",
                              "emo_reg_cluster_rMFG", "emo_reg_cluster_rPar", 
                              "emo_reg_cluster_rSFG",
                              "disgust_cluster_lAmyg", 
                              "disgust_cluster_lInsIFG", "disgust_cluster_lPar",
                              "disgust_cluster_rAmyg", "disgust_cluster_rClaus",
                              "disgust_cluster_rInsIFG",
                              "fear_cluster_lAmyg", 
                              "fear_cluster_lIFG", "fear_cluster_lThal", 
                              "fear_cluster_lmOccGyr", "fear_cluster_rAmyg", 
                              "fear_cluster_rIFG", "fear_cluster_rITG"),
                           labels=c("er_lCingG", "er_lIFG", "er_lMFG",
                                    "er_lMTG", "er_lSFG", "er_lTG", "er_rInsIFG",
                                    "er_rMFG", "er_rPar", "er_rSFG",
                                    "disgust_lAmyg", "disgust_lInsIFG", 
                                    "disgust_lPar", "disgust_rAmyg", "disgust_rClaus",
                                    "disgust_rInsIFG",
                                    "fear_lAmyg", "fear_lIFG", "fear_lThal", 
                                    "fear_lmOccGyr", "fear_rAmyg", "fear_rIFG", "fear_rITG"))

# Join with subject data
m3.use <- inner_join(x=m3.use, y=s.use, by = "participantID") 
#length(unique(m3.use$participantID)) #n = 42

# Clean the environment
rm(m1, m3)

#################################

#################################
# 2. Define Search Space Datasets ####
## Search Space ROIs ####
fspace <- c("fear_lAmyg", "fear_lIFG",
             "fear_lmOccGyr", "fear_lThal", "fear_rAmyg",
             "fear_rIFG", "fear_rITG")

dspace <- c("disgust_lAmyg", "disgust_lInsIFG",
             "disgust_lPar", "disgust_rAmyg", "disgust_rClaus", 
             "disgust_rInsIFG")

espace <- c("er_lCingG",
             "er_lIFG", "er_lMFG", "er_lMTG", 
             "er_lSFG", "er_lTG", "er_rInsIFG",
             "er_rMFG","er_rPar", "er_rSFG")

# naming differences for multivariate data
mv.dspace <- c("disgust_cluster_lAmyg", 
              "disgust_cluster_lInsIFG", "disgust_cluster_lPar",
              "disgust_cluster_rAmyg", "disgust_cluster_rClaus",
              "disgust_cluster_rInsIFG")

mv.fspace <- c("fear_cluster_lAmyg", 
                "fear_cluster_lIFG", "fear_cluster_lThal", 
                "fear_cluster_lmOccGyr", "fear_cluster_rAmyg", 
                "fear_cluster_rIFG", "fear_cluster_rITG") 

mv.espace <- c("emo_reg_cluster_lCingG", "emo_reg_cluster_lIFG",
                "emo_reg_cluster_lMFG", "emo_reg_cluster_lMTG", 
                "emo_reg_cluster_lSFG", "emo_reg_cluster_lTG",
                "emo_reg_cluster_rInsIFG", "emo_reg_cluster_rMFG",   
                "emo_reg_cluster_rPar", "emo_reg_cluster_rSFG")

## M1: Multivariate Fear & Disgust Spaces ####
# Reduce to ROIs needed to test M1
m1.wide.fd <- filter(m1.wide, roi %in% mv.dspace | roi %in% mv.fspace)

## M2: Univariate TOM Network (3 contrasts) ####
m2.tn <- m2 %>%
  filter(extraction_contrast == "con_1_tgn-gt-cgn")
m2.dn <- m2 %>%
  filter(extraction_contrast == "con_2_cgd-gt-cgn")
m2.fn <- m2 %>%
  filter(extraction_contrast == "con_3_cgf-gt-cgn")

## M3: Univariate Emotion Regulation ROI Dataset ####
m3.er <- filter(m3.use, roi %in% espace)
m3.er$roi_labels <- substring(m3.er$roi, first = 4)

## Manipulation Check: Univariate Fear ROI Dataset ####
mc.fear <- filter(m3.use, roi %in% fspace)
mc.fear$roi_labels <- substring(mc.fear$roi, first = 6)

## Manipulation Check: Univariate Disgust ROI Dataset ####
mc.disg <- filter(m3.use, roi %in% dspace)
mc.disg$roi_labels <- substring(mc.disg$roi, first = 9)


##################################

##################################
# 3. Normality Checks ####


## M1: Multivariate Data ####
### T Condition Correlations (Opponents) - Mixed ####
g1.tndf.shapiro <- m1.wide %>%
  filter(Group == 1) %>%
  select(participantID, roi, TD_r, TF_r, TN_r) %>%
  group_by(roi) %>%
  summarise(TN = list(TN_r),
            TD = list(TD_r),
            TF = list(TF_r)) %>% 
  group_by(roi) %>%
  mutate(TN_shapiro_p = shapiro.test(x = unlist(TN))$p.value,
         TD_shapiro_p = shapiro.test(x = unlist(TD))$p.value,
         TF_shapiro_p = shapiro.test(x = unlist(TF))$p.value) %>%
  select(-TN, -TD, -TF)
g1.tndf.shapiro

# TN/TD:
# - Normal pair: d_lAmyg, d_lInsIFG, d_lPar, d_rAmyg, d_rClaus, d_rInsIFG, 
# f_rAmyg, f_rIFG
# - Non-normal pair: f_lAmyg, f_lIFG, f_lmOccGyr, f_lThal, f_rITG.
# TN/TF: 
# - Normal pair: d_lAmyg, d_lInsIFG, d_lPar, d_rAmyg, d_rClaus, d_rInsIFG, 
# f_rAmyg, f_rIFG
# - Non-normal pair: f_lAmyg, f_lIFG, f_lmOccGyr, f_lThal, f_rITG
# TD/TF: 
# - Normal pair: d_lAmyg, d_lInsIFG, d_lPar, d_rAmyg, d_rClaus, d_rInsIFG, 
# f_lAmyg, f_rAmyg, f_rIFG
# - Non-normal pair:  f_lIFG, f_lmOccGyr, f_lThal, f_rITG

### T Condition correlations (Endorsers) - Mixed ####
g2.tndf.shapiro <- m1.wide %>%
  filter(Group == 2) %>%
  select(participantID, roi, TD_r, TF_r, TN_r) %>%
  group_by(roi) %>%
  summarise(TN = list(TN_r),
            TD = list(TD_r),
            TF = list(TF_r)) %>% 
  group_by(roi) %>%
  mutate(TN_shapiro_p = shapiro.test(x = unlist(TN))$p.value,
         TD_shapiro_p = shapiro.test(x = unlist(TD))$p.value,
         TF_shapiro_p = shapiro.test(x = unlist(TF))$p.value) %>%
  select(-TN, -TD, -TF)
g2.tndf.shapiro

# TN/TD:
# - Normal pair: d_lAmyg, d_lInsIFG, d_rAmyg, f_lAmyg
# f_lThal, f_rAmyg, f_rIFG
# - Non-normal pair: d_rClaus, d_rInsIFG, 
# f_lIFG, f_lmOccGyr,  f_rITG
# TN/TF: 
# - Normal pair: d_lAmyg, d_lInsIFG, d_lPar, d_rAmyg, d_rClaus, d_rInsIFG, 
# f_lAmyg,  f_lThal, f_rAmyg, f_rIFG
# - Non-normal pair: f_lIFG, f_lmOccGyr, f_rITG
# TD/TF: 
# - Normal pair: d_lAmyg, d_lInsIFG, d_lPar, d_rAmyg, d_rInsIFG, 
# f_lAmyg, f_lThal, f_rAmyg, f_rIFG,  
# - Non-normal pair: d_rClaus, 
# f_lIFG, f_lmOccGyr, f_rITG

### FDN Correlations - Mixed ####
fdn.grp.shapiro <- m1.wide.fd %>%
  select(participantID, Group, roi, wb_DF, wb_DN, wb_FN) %>%
  group_by(roi, Group) %>%
  summarise(wbDF = list(wb_DF),
            wbDN = list(wb_DN),
            wbFN = list(wb_FN)) %>% 
  group_by(roi, Group) %>%
  mutate(DF_shapiro_p = shapiro.test(x = unlist(wbDF))$p.value,
         DN_shapiro_p = shapiro.test(x = unlist(wbDN))$p.value,
         FN_shapiro_p = shapiro.test(x = unlist(wbFN))$p.value) %>%
  select(-wbDF, -wbDN, -wbFN)
fdn.grp.shapiro
# Opponents:
# DF/DN:
# - Normal pair:  f_lAmyg, f_lIFG, f_lThal
# - Non-normal pair: d_lAmyg, d_lInsIFG,  d_lPar, d_rAmyg, d_rClaus, 
# d_rInsIFG, f_lmOccGyr, f_rAmyg, f_rIFG, f_rITG
# DF/FN: 
# - Normal pair:   d_lInsIFG,  d_rInsIFG, f_lAmyg, f_lIFG, 
# - Non-normal pair: d_lAmyg, d_lPar, d_rAmyg, d_rClaus, 
# f_lmOccGyr, f_lThal, f_rAmyg, f_rIFG, f_rITG
# DN/FN: 
# - Normal pair: d_lInsIFG, f_lAmyg, f_lIFG  
# - Non-normal pair: d_lAmyg, d_lPar, d_lPar, d_rClaus, d_rInsIFG,
# f_lmOccGyr, f_lThal, f_rAmyg, f_rIFG, f_rITG

# Endorsers:
# DF/DN:
# - Normal pair:  d_lAmyg, f_lIFG, f_rAmyg
# - Non-normal pair: d_lInsIFG,  d_lPar, d_rAmyg, d_rClaus, 
# d_rInsIFG, f_lAmyg, f_lmOccGyr, f_lThal, f_rIFG, f_rITG
# DF/FN: 
# - Normal pair:   d_lInsIFG,  d_lPar, d_rInsIFG, 
# f_lAmyg, f_rIFG, f_rITG 
# - Non-normal pair: d_lAmyg, d_rAmyg, d_rClaus, f_lIFG, 
# f_lmOccGyr, f_lThal, f_rAmyg 
# DN/FN: 
# - Normal pair: f_lmOccGyr, f_lThal,  
# - Non-normal pair: d_lAmyg, d_lInsIFG,  d_lPar, d_rAmyg, d_rClaus, d_rInsIFG,
#  f_lAmyg, f_lIFG, f_rAmyg, f_rIFG, f_rITG

# Most pairs are non-normal; default testing is Yuen's D for dependent samples




## M2: ToM Network Data ####
### T>N Opponents - all normal ####
g1.tn.tom.shapiro <- m2.tn %>%
  filter(Group == 1) %>%
  group_by(ROI) %>%
  summarise(Mean_vox = list(mean_extracted_betas)) %>% 
  group_by(ROI) %>%
  mutate(TN_shapiro_p = round(shapiro.test(x = unlist(Mean_vox))$p.value, 4)) %>%
  select(-Mean_vox)
g1.tn.tom.shapiro
# all normal

### T>N Endorsers - all normal ####
g2.tn.tom.shapiro <- m2.tn %>%
  filter(Group == 2) %>%
  group_by(ROI) %>%
  summarise(Mean_vox = list(mean_extracted_betas)) %>% 
  group_by(ROI) %>%
  mutate(TN_shapiro_p = round(shapiro.test(x = unlist(Mean_vox))$p.value, 4)) %>%
  select(-Mean_vox)
g2.tn.tom.shapiro

### D>N Opponents - 2 not normal; dmPFC ok ####
g1.dn.tom.shapiro <- m2.dn %>%
  filter(Group == 1) %>%
  group_by(ROI) %>%
  summarise(Mean_vox = list(mean_extracted_betas)) %>% 
  group_by(ROI) %>%
  mutate(TN_shapiro_p = round(shapiro.test(x = unlist(Mean_vox))$p.value, 4)) %>%
  select(-Mean_vox)
g1.dn.tom.shapiro
# rTPJ and vmPFC not normal

### D>N Endorsers - 2 not normal, incl dmpfc ####
g2.dn.tom.shapiro <- m2.dn %>%
  filter(Group == 2) %>%
  group_by(ROI) %>%
  summarise(Mean_vox = list(mean_extracted_betas)) %>% 
  group_by(ROI) %>%
  mutate(TN_shapiro_p = round(shapiro.test(x = unlist(Mean_vox))$p.value, 4)) %>%
  select(-Mean_vox)
g2.dn.tom.shapiro
# dmPFC and vmPFC not normal

### F>N Opponents - all normal ####
g1.fn.tom.shapiro <- m2.fn %>%
  filter(Group == 1) %>%
  group_by(ROI) %>%
  summarise(Mean_vox = list(mean_extracted_betas)) %>% 
  group_by(ROI) %>%
  mutate(TN_shapiro_p = round(shapiro.test(x = unlist(Mean_vox))$p.value, 4)) %>%
  select(-Mean_vox)
g1.fn.tom.shapiro
# all normal

### F>N Endorsers - all normal ####
g2.fn.tom.shapiro <- m2.fn %>%
  filter(Group == 2) %>%
  group_by(ROI) %>%
  summarise(Mean_vox = list(mean_extracted_betas)) %>% 
  group_by(ROI) %>%
  mutate(TN_shapiro_p = round(shapiro.test(x = unlist(Mean_vox))$p.value, 4)) %>%
  select(-Mean_vox)
g2.fn.tom.shapiro
# all normal



## M3 + Other Tests: CBMA Univariate Contrasts ####
### ER Space Only: T > N (Opponents) - 3 non-normal ####
g1.tn.univ.shapiro <- m3.use %>%
  filter(Group == 1) %>%
  filter(roi %in% espace, contrast == "tgn-gt-cgn") %>%
  group_by(roi) %>%
  summarise(Mean_vox = list(mean_topvoxels)) %>% 
  group_by(roi) %>%
  mutate(TN_shapiro_p = round(shapiro.test(x = unlist(Mean_vox))$p.value, 4)) %>%
  select(-Mean_vox)
g1.tn.univ.shapiro
# Non-normal: er_lCing, er_rInsIFG, er_lTG

### ER Space Only: T > N (Endorsers) - 1 non-normal ####
g2.tn.univ.shapiro <- m3.use %>%
  filter(Group == 2) %>%
  filter(roi %in% espace, contrast == "tgn-gt-cgn") %>%
  group_by(roi) %>%
  summarise(Mean_vox = list(mean_topvoxels)) %>% 
  group_by(roi) %>%
  mutate(TN_shapiro_p = round(shapiro.test(x = unlist(Mean_vox))$p.value, 4)) %>%
  select(-Mean_vox)
g2.tn.univ.shapiro
# Non-normal: er_rInsIFG


### All CBMA Spaces: F > N (Opponents) - 2 not normal ####
g1.fn.univ.shapiro <- m3.use %>%
  filter(Group == 1) %>%
  filter(contrast == "cgf-gt-cgn") %>%
  group_by(roi) %>%
  summarise(Mean_vox = list(mean_topvoxels)) %>% 
  group_by(roi) %>%
  mutate(FN_shapiro_p = round(shapiro.test(x = unlist(Mean_vox))$p.value, 4)) %>%
  select(-Mean_vox)
g1.fn.univ.shapiro

### All CBMA Spaces: F > N (Endorsers) - 1 non-normal ####
g2.fn.univ.shapiro <- m3.use %>%
  filter(Group == 2) %>%
  filter(contrast == "cgf-gt-cgn") %>%
  group_by(roi) %>%
  summarise(Mean_vox = list(mean_topvoxels)) %>% 
  group_by(roi) %>%
  mutate(FN_shapiro_p = round(shapiro.test(x = unlist(Mean_vox))$p.value, 4)) %>%
  select(-Mean_vox)
g2.fn.univ.shapiro
# All except er_lMTG are normal

### All CBMA Spaces: F > N (all subjects) - 5 non-normal ####
fn.univ.shapiro <- m3.use %>%
  filter(contrast == "cgf-gt-cgn") %>%
  group_by(roi) %>%
  summarise(Mean_vox = list(mean_topvoxels)) %>% 
  group_by(roi) %>%
  mutate(FN_shapiro_p = round(shapiro.test(x = unlist(Mean_vox))$p.value, 4)) %>%
  select(-Mean_vox)
fn.univ.shapiro

# Adjust for fear_lmOccGyr in hm1a test
### All CBMA Spaces: D > N (Opponents) - 2 non-normal ####
g1.dn.univ.shapiro <- m3.use %>%
  filter(Group == 1) %>%
  filter(contrast == "cgd-gt-cgn") %>%
  group_by(roi) %>%
  summarise(Mean_vox = list(mean_topvoxels)) %>% 
  group_by(roi) %>%
  mutate(DN_shapiro_p = round(shapiro.test(x = unlist(Mean_vox))$p.value, 4)) %>%
  select(-Mean_vox)
g1.dn.univ.shapiro
# Everything is normal except er_lMFG, er_lMTG

### All CBMA Spaces: D > N (Endorsers) - 1 non-normal ####
g2.dn.univ.shapiro <- m3.use %>%
  filter(Group == 2) %>%
  filter(contrast == "cgd-gt-cgn") %>%
  group_by(roi) %>%
  summarise(Mean_vox = list(mean_topvoxels)) %>% 
  group_by(roi) %>%
  mutate(DN_shapiro_p = round(shapiro.test(x = unlist(Mean_vox))$p.value, 4)) %>%
  select(-Mean_vox)
g2.dn.univ.shapiro
# Everything is normal except er_rPar

### All CBMA Spaces: D > N (all subjects) - 2 non-normal ####
dn.univ.shapiro <- m3.use %>%
  filter(contrast == "cgd-gt-cgn") %>%
  group_by(roi) %>%
  summarise(Mean_vox = list(mean_topvoxels)) %>% 
  group_by(roi) %>%
  mutate(DN_shapiro_p = round(shapiro.test(x = unlist(Mean_vox))$p.value, 4)) %>%
  select(-Mean_vox)
dn.univ.shapiro

# No adjustments needed for HM2a test

##################################

##################################
# 4. M1 (Threat Emotions) Tests ####
# Note that corr comparisons are dependent, so
# Yuen's D is used (t-test for dependent non-normal w/ 20% default mean trim); 
# bootstrap is used for p-values (WRS function is ydbt(), seed is set by default)
# When groups are independent, yuenbt() is used
# When 1-sample, 1-sided tests are required, the custom function trdiff1s is used
# (same as the 1-sample trimmed mean function trimcibt, except p-values are 
# only calculated for the case X > 0, i.e., alternative = "g")

# Steps as discussed in pre-registration

## Step 1: Identify ROIs with Double Dissociation (F vs. D vs. N)  ####

fdn.diss <- m1.wide.fd %>%
  select(participantID, roi, wb_DF, wb_DN, wb_FN) %>%
  group_by(roi) %>%
  summarise(wbDF = list(wb_DF),
            wbDN = list(wb_DN),
            wbFN = list(wb_FN)) %>% 
  group_by(roi) %>%
  mutate(FD_trm_tr1s = trdiff1sp(x = unlist(wbDF))$estimate,
         FD_stat_tr1s = trdiff1sp(x = unlist(wbDF))$test.stat,
         FD_pvalue_tr1s = trdiff1sp(x = unlist(wbDF))$p.value,
         FD_n_tr1s = trdiff1sp(x = unlist(wbDF))$n,
         DN_trm_tr1s = trdiff1sp(x = unlist(wbDN))$estimate,
         DN_stat_tr1s = trdiff1sp(x = unlist(wbDN))$test.stat,
         DN_pvalue_tr1s = trdiff1sp(x = unlist(wbDN))$p.value,
         DN_n_tr1s = trdiff1sp(x = unlist(wbDN))$n,
         FN_trm_tr1s = trdiff1sp(x = unlist(wbFN))$estimate,
         FN_stat_tr1s = trdiff1sp(x = unlist(wbFN))$test.stat,
         FN_pvalue_tr1s = trdiff1sp(x = unlist(wbFN))$p.value,
         FN_n_tr1s = trdiff1sp(x = unlist(wbFN))$n) %>%
  select(-wbDF, -wbDN, -wbFN)
fdn.diss
# This df is the input for tables C6 & C7 in Appendix C

# Double-Dissociations:
# Dissociation p-value (not adjusted for number of tests)
dissp <- 0.05

f.test <- ifelse(fdn.diss$FD_pvalue_tr1s < dissp &
                        fdn.diss$FN_pvalue_tr1s < dissp, 1, 0)
d.test <- ifelse(fdn.diss$FD_pvalue_tr1s < dissp &
                      fdn.diss$DN_pvalue_tr1s < dissp, 1, 0)
fdn.use <- data.frame(roi = fdn.diss$roi,
                        f_diss = f.test,
                        d_diss = d.test)
fdn.use

# Define reduced search spaces as ROIs with double dissociations
m1.f.ss <- filter(fdn.use, f_diss == 1 & roi %in% mv.fspace) %>%
  select(roi) %>% c() %>% unlist() %>% as.character()
m1.f.ss
# 3 fear w 1000 bootstrap samples 
# "fear_cluster_lmOccGyr" "fear_cluster_rAmyg" "fear_cluster_rITG"    

m1.d.ss <- filter(fdn.use, d_diss == 1 & roi %in% mv.dspace) %>%
  select(roi) %>% c() %>% unlist() %>% as.character()
m1.d.ss
# All six regions pass double dissociation
# "disgust_cluster_lAmyg"   "disgust_cluster_lInsIFG"
# "disgust_cluster_lPar"    "disgust_cluster_rAmyg"  
# "disgust_cluster_rClaus"  "disgust_cluster_rInsIFG"

# Combine into M1 Testing search space
diss.rois <- c(m1.f.ss, m1.d.ss)

## Step 2: Significance Testing in Endorsers (H1a Tests) ####
# Per normality checks for Endorsers: 
# TN/TD:
# - Normal pair: d_lAmyg, d_lInsIFG, d_lPar, f_rAmyg
# - Non-normal pair: f_rITG, d_rInsIFG
# TN/TF: 
# - Normal pair: d_lAmyg, d_lInsIFG, d_lPar, d_rInsIFG, 
# f_rAmyg
# - Non-normal pair: f_rITG
# TD/TF: 
# - Normal pair: d_lAmyg, d_lInsIFG, d_lPar, d_rInsIFG,
# f_rAmyg  
# - Non-normal pair: f_rITG

# Yuen's Robust T test for Dependent Groups (trim = 0.2)
# Note: Most ROIs are normal, but paper presents robust results to be conservative
# Take-aways don't change for non-robust version of test (trim = 0) 
h1a.trim = 0.2

set.seed(1234)
h1a.m1.ydbt <- m1.wide.fd %>%
  filter(Group == 2) %>%
  # search only in confirmed double dissociation ROIS
  filter(roi %in% diss.rois) %>%
  select(participantID, roi, TD_r, TF_r, TN_r) %>%
  group_by(roi) %>%
  summarise(TD = list(TD_r),
            TF = list(TF_r),
            TN = list(TN_r)) %>%
  group_by(roi) %>%
    mutate(TD_trm = mean(unlist(TD), trim = h1a.trim),
           TF_trm = mean(unlist(TF), trim = h1a.trim),
           TN_trm = mean(unlist(TN), trim = h1a.trim),
           TNTD_diff = yuend(x = unlist(TD), y = unlist(TN), tr = h1a.trim)$diff,
           TNTD_pvalue = yuend(x = unlist(TD), y = unlist(TN), tr = h1a.trim, nboot = 1000)$p.value,
           TNTF_diff = yuend(x = unlist(TF), y = unlist(TN), tr = h1a.trim)$diff,
           TNTF_pvalue = yuend(x = unlist(TN), y = unlist(TF), tr = h1a.trim, nboot = 1000)$p.value,
           TDTF_diff = yuend(x = unlist(TD), y = unlist(TF), tr = h1a.trim)$diff,
           TDTF_pvalue = yuend(x = unlist(TD), y = unlist(TF), tr = h1a.trim, nboot = 1000)$p.value) %>%
  select(-TD, -TF, -TN)
print(h1a.m1.ydbt)

# Test for any ROI where T is significantly closer to 
# D (or F) than to F (or D) or N 
# Such that corr w target emotion (TD or TF)
# is greater than TD/TF and TN

# H1a expectation is T is closer to D than F or N
h1a.m1.disgust <- select(h1a.m1.ydbt, roi,
                     TNTD_diff, TDTF_diff,
                     TNTD_pvalue, TDTF_pvalue) %>% 
  filter(roi %in% mv.dspace)
# Reverse TNTD_diff for ease of interpretation
h1a.m1.disgust$TDTN_diff <- -h1a.m1.disgust$TNTD_diff
h1a.m1.disgust <- select(h1a.m1.disgust, -TNTD_diff)
# Run tests for significant proximity
h1a.m1.disgust$DvsN <- ifelse(h1a.m1.disgust$TDTN_diff > 0 &
                                h1a.m1.disgust$TNTD_pvalue < 0.05, 1, 0)
h1a.m1.disgust$DvsF <- ifelse(h1a.m1.disgust$TDTF_diff > 0 &
                                h1a.m1.disgust$TDTF_pvalue < 0.05, 1, 0)
h1a.m1.disgust$D_sim <- ifelse(h1a.m1.disgust$DvsN == 1 & 
                                 h1a.m1.disgust$DvsF == 1, 1, 0)
h1a.m1.disgust
# No regions pass test of double dissociation

# Not included, but checked 
h1a.m1.fear <- select(h1a.m1.ydbt, roi,
                      TNTF_diff, TDTF_diff,
                      TNTF_pvalue, TDTF_pvalue) %>% 
  filter(roi %in% mv.fspace)
# Reverse TNTF_diff and TDTF_diff for ease of interpretation
h1a.m1.fear$TFTN_diff <- -h1a.m1.fear$TNTF_diff
h1a.m1.fear$TFTD_diff <- -h1a.m1.fear$TDTF_diff
h1a.m1.fear <- select(h1a.m1.fear, -TNTF_diff, -TDTF_diff)
# Run tests for significant proximity
h1a.m1.fear$FvsN <- ifelse(h1a.m1.fear$TFTN_diff > 0 &
                                h1a.m1.fear$TNTF_pvalue < 0.05, 1, 0)
h1a.m1.fear$FvsD <- ifelse(h1a.m1.fear$TFTD_diff > 0 &
                                h1a.m1.fear$TDTF_pvalue < 0.05, 1, 0)
h1a.m1.fear$F_sim <- ifelse(h1a.m1.fear$FvsN == 1 & 
                                 h1a.m1.fear$FvsD == 1, 1, 0)
h1a.m1.fear
# No regions pass test of double dissociation


# NOTE: These are the inputs for Tables D1-D2 in Appendix D


## Opponents (H1b Tests) ####

# Analyses are not limited to regions of dissociation 
# About 50% of ROIs are non-normal per tests above and paper
# reports robust version only (trim = 0.2)
# Take-aways unchanged for regions where normality is assumed (trim = 0)
g1fd.trim = 0.2

set.seed(1234)
g1.fd.m1.ydbt <- m1.wide.fd %>%
  filter(Group == 1) %>%
  select(participantID, roi, TD_r, TF_r, TN_r) %>%
  group_by(roi) %>%
  summarise(TD = list(TD_r),
            TF = list(TF_r),
            TN = list(TN_r)) %>% 
  group_by(roi) %>%
  mutate(TD_trm = mean(unlist(TD), trim = g1fd.trim),
         TF_trm = mean(unlist(TF), trim = g1fd.trim),
         TN_trm = mean(unlist(TN), trim = g1fd.trim),
         TNTD_diff = yuend(x = unlist(TN), y = unlist(TD), tr = g1fd.trim)$diff,
         TNTD_pvalue = yuend(x = unlist(TN), y = unlist(TD), tr = g1fd.trim, nboot = 1000)$p.value,
         TNTF_diff = yuend(x = unlist(TN), y = unlist(TF), tr = g1fd.trim)$diff,
         TNTF_pvalue = yuend(x = unlist(TN), y = unlist(TF), tr = g1fd.trim, nboot = 1000)$p.value,
         TDTF_diff = yuend(x = unlist(TD), y = unlist(TF), tr = g1fd.trim)$diff,
         TDTF_pvalue = yuend(x = unlist(TD), y = unlist(TF), tr = g1fd.trim, nboot = 1000)$p.value) %>%
  select(-TD, -TF, -TN)
print(g1.fd.m1.ydbt)

# Test for any ROI where T is significantly closer to 
# N than to D or F
# In disgust ROIs
h1b.m1.disgust <- select(g1.fd.m1.ydbt, roi,
                         TNTD_diff, TNTF_diff,
                         TNTD_pvalue, TNTF_pvalue) %>% 
  filter(roi %in% mv.dspace)
# Run tests for significant proximity
h1b.m1.disgust$NvsD <- ifelse(h1b.m1.disgust$TNTD_diff > 0 &
                                h1b.m1.disgust$TNTD_pvalue < 0.05, 1, 0)
h1b.m1.disgust$NvsF <- ifelse(h1b.m1.disgust$TNTF_diff > 0 &
                                h1b.m1.disgust$TNTF_pvalue < 0.05, 1, 0)
h1b.m1.disgust$N_sim <- ifelse(h1b.m1.disgust$NvsD == 1 & 
                                 h1b.m1.disgust$NvsF == 1, 1, 0)
h1b.m1.disgust
# No regions pass test of double dissociation

# In fear ROIs
h1b.m1.fear <- select(g1.fd.m1.ydbt, roi,
                      TNTD_diff, TNTF_diff,
                      TNTD_pvalue, TNTF_pvalue) %>% 
  filter(roi %in% mv.fspace)
# Run tests for significant proximity
h1b.m1.fear$NvsD <- ifelse(h1b.m1.fear$TNTD_diff > 0 &
                                h1b.m1.fear$TNTD_pvalue < 0.05, 1, 0)
h1b.m1.fear$NvsF <- ifelse(h1b.m1.fear$TNTF_diff > 0 &
                                h1b.m1.fear$TNTF_pvalue < 0.05, 1, 0)
h1b.m1.fear$N_sim <- ifelse(h1b.m1.fear$NvsD == 1 & 
                                 h1b.m1.fear$NvsF == 1, 1, 0)
h1b.m1.fear
# No regions pass test of double dissociation

# This is the input for Tables 3-4 in Appendix D


## Group Differences (H1c Tests) ####
g1_tn_fd <- m1.wide.fd %>%
  filter(Group == 1) %>%
  select(participantID, roi, TN_r) %>%
  group_by(roi) %>%
  summarise(NT_g1 = list(TN_r))

g2_tn_fd <- m1.wide.fd %>%
  filter(Group == 2) %>%
  select(participantID, roi, TN_r) %>%
  group_by(roi) %>%
  summarise(NT_g2 = list(TN_r))

all_tn_fd <- inner_join(g1_tn_fd, g2_tn_fd, by = "roi")

# Mixed results for normality so we use Yuen's robust t-test 
# for independent samples (trim = 0.2)

# list to hold output
h1c.m1.ls <- list()
# trim to use for robust version
h1c.trim = 0.2
# 0.2 is default in WRS2 pkg

set.seed(1234)
for(i in 1:length(all_tn_fd$roi)){
  # extract the row
  roi.use <- all_tn_fd[i,]
  # unlist the variables 
  g1.use <- unlist(roi.use$NT_g1)
  g2.use <- unlist(roi.use$NT_g2)
  # combine into dataframe
  roi.df <- data.frame(NT_g1 = g1.use,
                       NT_g2 = g2.use)
  # yuen takes a long df
  roi.df.long <- pivot_longer(roi.df, cols = everything(), names_to = "Group", values_to = "NT_sim")
  roi.ytest <- yuenbt(NT_sim ~ Group, data = roi.df.long, tr = h1c.trim, nboot = 1000)
  # trimmed means - just to check same as ytest
  g1_mean <- mean(roi.df$NT_g1, trim=h1c.trim)
  g2_mean <- mean(roi.df$NT_g2, trim=h1c.trim)
  # Save outputs 
  h1c.m1.ls[[i]] <- data.frame(roi = roi.use$roi,
                               g1_tm = mean(g1.use, trim = h1c.trim),
                               g2_tm = mean(g2.use, trim = h1c.trim),
                               y_diff = roi.ytest$diff,
                       y_stat = unname(roi.ytest$test),
                       y_pval = unname(roi.ytest$p.value))
}

h1c.m1.df <- bind_rows(h1c.m1.ls)
rownames(h1c.m1.df) <- NULL 
h1c.m1.df

# No regions of difference between groups

###################################

###################################
# 5. M2 (Mentalizing) Tests ####
### Opponents (H2a Tests) ####
# All data assumed normal (see above)
h2a.test <- m2.tn %>%
  filter(Group == 1) %>%
  group_by(ROI) %>%
  group_map(~ t.test(.$mean_extracted_betas,
                     mu = 0,
                     alternative = "greater"))
names(h2a.test) <- unique(m2.tn$ROI)

h2a.test.df <- lapply(h2a.test, function(x) x[c('estimate', 'statistic', 'p.value')])
h2a.test.df <- bind_rows(h2a.test.df)
h2a.test.df$ROI <- c("DMPFC", "l. TPJ",
                     "MMPFC", "Precuneus",
                     "r. STS", "r. TPJ",
                     "VMPFC")
h2a.test.df
# No evidence of network-wide activation

# ROI results reported in MS
h2a.test[["dmPFC"]] # ***
#M = 0.9105438; t = 5.145, df = 20, p-value = 2.465e-05

h2a.test[["mmPFC"]] # **
#M = 0.6476285; t = 2.9578, df = 20, p-value = 0.00389


### Endorsers (same test for reference) ####
# No direction specified; 
h2a.check <- m2.tn %>%
  filter(Group == 2) %>%
  group_by(ROI) %>%
  group_map(~ t.test(.$mean_extracted_betas,
                     mu = 0,
                     alternative = "two.sided"))
names(h2a.check) <- unique(m2.tn$ROI)

h2a.check.df <- lapply(h2a.check, function(x) x[c('estimate', 'statistic', 'p.value')])
h2a.check.df <- bind_rows(h2a.check.df)
h2a.check.df$ROI <- c("DMPFC", "l. TPJ",
                     "MMPFC", "Precuneus",
                     "r. STS", "r. TPJ",
                     "VMPFC")
h2a.check.df
# No evidence of network-wide activation

# Reported in paper:
h2a.check[["dmPFC"]] 
#M = 0.1082918; t = 0.54417, df = 19, p-value = 0.5927

h2a.check[["mmPFC"]]
#M = 0.002547675; t = 0.009061, df = 19, p-value = 0.9929


### Between-Groups (H2b Tests) ####
# All SW tests were normal for both groups
# so Welch's T-Test used (two-sided, per H2b which was for difference not direction)
# Group 1 is always first in test;

# Run a 2sample t-test in each ROI
h2b.tests <- grpttest.tom(x = m2.tn,
                          z = "two.sided")
names(h2b.tests) <- unique(m2.tn$ROI)

# Summarize 
h2b.means.df <- lapply(h2b.tests, function(x) x[c('estimate')])
h2b.mean.g1.df <- lapply(h2b.means.df, function(x) x$estimate[1])
h2b.mean.g1.df <- bind_rows(h2b.mean.g1.df)
h2b.mean.g2.df <- lapply(h2b.means.df, function(x) x$estimate[2])
h2b.mean.g2.df <- bind_rows(h2b.mean.g2.df)

h2b.stats.df <- lapply(h2b.tests, function(x) x[c('statistic', 'p.value')])
h2b.stats.df <- bind_rows(h2b.stats.df)

h2b.tests.df <- data.frame(mean_g1 = h2b.mean.g1.df,
                           mean_g2 = h2b.mean.g2.df,
                           t_statistic = h2b.stats.df[,1],
                           p_value = h2b.stats.df[,2])
h2b.tests.df$ROI <- c("DMPFC", "l. TPJ",
                      "MMPFC", "Precuneus",
                      "r. STS", "r. TPJ",
                      "VMPFC")
h2b.tests.df

# Reported in paper:
h2b.tests[["dmPFC"]] # **
#t = 3.0124, df = 38.224, p-value = 0.004579

h2b.tests[["mmPFC"]]
#t = 1.8101, df = 36.337, p-value = 0.07855


### DMPFC Negative Emotion Check ####

# Disgust 
# Need to use robust method to check on dmPFC per SW testing
m2.dmpfc.dn <- filter(m2.dn, ROI == "dmPFC")
str(m2.dmpfc.dn)

# Without bootstrap
YuenTTest(filter(m2.dmpfc.dn, Group == 1)$mean_extracted_betas,
          filter(m2.dmpfc.dn, Group == 2)$mean_extracted_betas,
          alternative = "two.sided")
#t = 0.6264, df = 22.761, trim = 0.200, p-value = 0.5373

# With bootstrap - reported
set.seed(1234)
yuenbt(mean_extracted_betas ~ Group, data = m2.dmpfc.dn, nboot = 1000, side = T)
# Test statistic: 0.6465 (df = NA), p-value = 0.521

# Fear
# Both groups assumed normal per SW testing
m2.dmpfc.fn <- filter(m2.fn, ROI == "dmPFC")

t.test(filter(m2.dmpfc.fn, Group == 1)$mean_extracted_betas,
       filter(m2.dmpfc.fn, Group == 2)$mean_extracted_betas,
       alternative = "two.sided",
       paired = F)
# t = -1.0644, df = 34.271, p-value = 0.2946

### Repeat with matched N=40 ####

# Drop the sub whose match failed ToM
ids42 <- unique(m1.wide$participantID)
ids41 <- unique(m2.tn$participantID)

id.missing <- setdiff(ids42, ids41)
# missing = sub-SAXEEMOfd29

# Sub 29's match is sub 35
id.drop <- c("sub-SAXEEMOfd35")

m2.tn.40 <-filter(m2.tn, !participantID == id.drop)

h2a.test.n40 <- m2.tn.40 %>%
  filter(Group == 1) %>%
  group_by(ROI) %>%
  group_map(~ t.test(.$mean_extracted_betas,
                     mu = 0,
                     alternative = "greater"))
names(h2a.test.n40) <- unique(m2.tn$ROI)

h2a.test.n40.df <- lapply(h2a.test.n40, function(x) x[c('estimate', 'statistic', 'p.value')])
h2a.test.n40.df <- bind_rows(h2a.test.n40.df)
h2a.test.n40.df$ROI <- c("DMPFC", "l. TPJ",
                         "MMPFC", "Precuneus",
                         "r. STS", "r. TPJ",
                         "VMPFC")
h2a.test.n40.df
# Same substantive results as N=41

h2b.tests.n40 <- grpttest.tom(x = m2.tn.40,
                              z = "two.sided")
names(h2b.tests.n40) <- unique(m2.tn$ROI)

# Summarize 
h2b.means.n40.df <- lapply(h2b.tests.n40, function(x) x[c('estimate')])
h2b.mean.n40.g1.df <- lapply(h2b.means.n40.df, function(x) x$estimate[1])
h2b.mean.n40.g1.df <- bind_rows(h2b.mean.n40.g1.df)
h2b.mean.n40.g2.df <- lapply(h2b.means.n40.df, function(x) x$estimate[2])
h2b.mean.n40.g2.df <- bind_rows(h2b.mean.n40.g2.df)

h2b.stats.n40.df <- lapply(h2b.tests.n40, function(x) x[c('statistic', 'p.value')])
h2b.stats.n40.df <- bind_rows(h2b.stats.n40.df)

h2b.tests.n40.df <- data.frame(mean_g1 = h2b.mean.n40.g1.df,
                               mean_g2 = h2b.mean.n40.g2.df,
                               t_statistic = h2b.stats.n40.df[,1],
                               p_value = h2b.stats.n40.df[,2])
h2b.tests.n40.df$ROI <- c("DMPFC", "l. TPJ",
                          "MMPFC", "Precuneus",
                          "r. STS", "r. TPJ",
                          "VMPFC")
h2b.tests.n40.df
# Same substantive results as N=41



###################################

###################################
# 6. M3 Testing: Self-Regulation ####

# Based on normality testing, robust t-tests for: 
# - T>N, Opponents: er_lCing, er_rInsIFG, er_lTG
# - T>N, Endorsers: er_rInsIFG

## Opponents (H3a Tests) ####

### Robust One Sample T-Test - Reported ####
# Required for er_lCingG, er_lTG, er_rInsIFG, per normality checks
# Reporting robust as it is more conservative, but 
# Standard t-test valid for 7 regions (see below; all significance tests yield same conclusion)

h3a.test <- filter(m3.er, 
                   Group == 1,
                   contrast == "tgn-gt-cgn") %>%
  group_by(roi) %>%
  group_map(~ trdiff1sp(.$mean_topvoxels))
names(h3a.test) <- unique(m3.er$roi)

# Convert results to table
h3a.test.est <- lapply(h3a.test, function(x) x$estimate)
h3a.test.stats <- lapply(h3a.test, function(x) x[c('test.stat', 'p.value')])

h3a.test.est.df <- unlist(h3a.test.est)
h3a.test.stats.df <- bind_rows(h3a.test.stats)

h3a.test.df <- data.frame(estimate = h3a.test.est.df,
                          t.stat = h3a.test.stats.df$test.stat,
                          p.value = h3a.test.stats.df$p.value)
h3a.test.df$ROI <- c("l. Cingulate", "l. IFG", "l. MFG", 
                     "l. MTG", "l. SFG", "l. TG", "rInsIFG", 
                     "r. MFG", "r. Parietal","r. SFG")
row.names(h3a.test.df) <- NULL
h3a.test.df


### Standard T-Tests (7 ROI) - not reported ####
h3a.test.std <- m3.er %>%
  filter(contrast == "tgn-gt-cgn", Group == 1) %>%
  group_by(roi) %>%
  group_map(~ t.test(.$mean_topvoxels,
                     mu = 0,
                     alternative = "greater"))
names(h3a.test.std) <- unique(m3.er$roi)


# Convert results to table
h3a.test.std.est <- lapply(h3a.test.std, function(x) x$estimate)
h3a.test.std.stats <- lapply(h3a.test.std, function(x) x[c('statistic', 'p.value')])

h3a.test.std.est.df <- unlist(h3a.test.std.est)
h3a.test.std.stats.df <- bind_rows(h3a.test.std.stats)

h3a.test.std.df <- data.frame(estimate = h3a.test.std.est.df,
                          t.stat = h3a.test.std.stats.df$statistic,
                          p.value = h3a.test.std.stats.df$p.value)

# drop er_lCingG, er_lTG, er_rInsIFG (requires robust)
h3a.test.std.df <- h3a.test.std.df[-c(1, 6, 7),]

h3a.test.std.df$ROI <- c("l. IFG", "l. MFG", 
                     "l. MTG", "l. SFG", 
                     "r. MFG", "r. Parietal","r. SFG")
h3a.test.std.df
# All significance tests hold  

## Endorsers (same tests for reference) ####

### Robust One Sample T-Test - Reported ####
# Only required for er_rInsIFG, but reported for comparability
# Standard t-tests replicate the same null results in 9 valid ROIS (see below)

h3a.check <- filter(m3.er, 
                    Group == 2,
                    contrast == "tgn-gt-cgn") %>%
  group_by(roi) %>%
  group_map(~ trdiff1sp(.$mean_topvoxels))
names(h3a.check) <- unique(m3.er$roi)

h3a.check.est <- lapply(h3a.check, function(x) x$estimate)
h3a.check.stats <- lapply(h3a.check, function(x) x[c('test.stat', 'p.value')])

h3a.check.est.df <- unlist(h3a.check.est)
h3a.check.stats.df <- bind_rows(h3a.check.stats)

h3a.check.df <- data.frame(estimate = h3a.check.est.df,
                          t.stat = h3a.check.stats.df$test.stat,
                          p.value = h3a.check.stats.df$p.value)
h3a.check.df$ROI <- c("l. Cingulate", "l. IFG", "l. MFG", 
                     "l. MTG", "l. SFG", "l. TG", "rInsIFG", 
                     "r. MFG", "r. Parietal","r. SFG")
row.names(h3a.check.df) <- NULL
h3a.check.df

### Standard T-Tests (9 ROI) - not reported ####
h3a.check.std <- m3.er %>%
  filter(contrast == "tgn-gt-cgn", Group == 2) %>%
  group_by(roi) %>%
  group_map(~ t.test(.$mean_topvoxels,
                     mu = 0,
                     alternative = "greater"))
names(h3a.check.std) <- unique(m3.er$roi)

# Convert results to table
h3a.check.std.est <- lapply(h3a.check.std, function(x) x$estimate)
h3a.check.std.stats <- lapply(h3a.check.std, function(x) x[c('statistic', 'p.value')])

h3a.check.std.est.df <- unlist(h3a.check.std.est)
h3a.check.std.stats.df <- bind_rows(h3a.check.std.stats)

h3a.check.std.df <- data.frame(estimate = h3a.check.std.est.df,
                              t.stat = h3a.check.std.stats.df$statistic,
                              p.value = h3a.check.std.stats.df$p.value)

# drop er_rInsIFG (requires robust)
h3a.check.std.df <- h3a.check.std.df[-c(7),]

h3a.check.std.df$ROI <- c(". Cingulate", "l. IFG", "l. MFG", 
                         "l. MTG", "l. SFG", "l. TG", 
                         "r. MFG", "r. Parietal","r. SFG")
h3a.check.std.df
# All significance tests hold (null results)




## Between-Groups (H3b Tests) ####

### Robust Two-Sample Tests - Reported ####
# Required for er_lCingG, er_lTG, er_rInsIFG
# Make sure alternative = "greater" in helper function

h3b.test <- grp.robust(m3.er, "tgn-gt-cgn")
names(h3b.test) <- unique(m3.er$roi)

# Convert results to table
h3b.test.mean1 <- lapply(h3b.test, function(x) x$estimate[1])
h3b.test.mean2 <- lapply(h3b.test, function(x) x$estimate[2])
h3b.test.stats <- lapply(h3b.test, function(x) x[c('statistic', 'p.value')])

h3b.test.mean1.df <- unlist(h3b.test.mean1)
h3b.test.mean2.df <- unlist(h3b.test.mean2)
h3b.test.stats.df <- bind_rows(h3b.test.stats)

h3b.test.df <- data.frame(g1.mean = h3b.test.mean1.df,
                          g2.mean = h3b.test.mean2.df,
                          t.stat = h3b.test.stats.df$statistic,
                          p.value = h3b.test.stats.df$p.value)
h3b.test.df$ROI <- c("l. Cingulate", "l. IFG", "l. MFG", 
                     "l. MTG", "l. SFG", "l. TG", "rInsIFG", 
                     "r. MFG", "r. Parietal","r. SFG")
row.names(h3b.test.df) <- NULL
h3b.test.df




### Standard 2-Sample T-tests (7 ROI) - reported for l. IFG only ####
# Welch's 2 Sample T-Tests inside custom function
h3b.test.std <- grpttest(m3.er, 
                       "tgn-gt-cgn", 
                       "greater")
names(h3b.test.std) <- unique(m3.er$roi)

h3b.test.std.mean1 <- lapply(h3b.test.std, function(x) x$estimate[1])
h3b.test.std.mean2 <- lapply(h3b.test.std, function(x) x$estimate[2])
h3b.test.std.stats <- lapply(h3b.test.std, function(x) x[c('statistic', 'p.value')])

h3b.test.std.mean1.df <- unlist(h3b.test.std.mean1)
h3b.test.std.mean2.df <- unlist(h3b.test.std.mean2)
h3b.test.std.stats.df <- bind_rows(h3b.test.std.stats)

h3b.test.std.df <- data.frame(g1.mean = h3b.test.std.mean1.df,
                              g2.mean = h3b.test.std.mean2.df,
                              t.stat = h3b.test.std.stats.df$statistic,
                              p.value = h3b.test.std.stats.df$p.value)

# drop er_lCingG, er_lTG, er_rInsIFG (requires robust)
h3b.test.std.df <- h3b.test.std.df[-c(1, 6, 7),]
h3b.test.std.df$ROI <- c("l. IFG", "l. MFG", 
                         "l. MTG", "l. SFG", 
                         "r. MFG", "r. Parietal","r. SFG")
h3b.test.std.df
# All significance tests hold  



############################################

############################################
# 7. Additional Pre-Registered Hypothesis tests ####
# T-test is the pre-registered method, but for   
# l. med. occipital gyrus (non-normal) we report the  
# robust method (trimmed mean, bootstrapped p-values)
# Note that whole-brain results (H-M1b, H-M2b) are reported in tables E1 and E2

## Fear Manipulation (H-M1a) #### 
### F > N (All Subjects), One Sample T-Test - Reported for 6 ROI ####
hm1a.std.roi <- c("fear_lAmyg", "fear_lIFG" ,"fear_lThal", 
                  "fear_rAmyg", "fear_rIFG", "fear_rITG")                       
  
hm1a.test.std <- mc.fear %>%
  filter(contrast == "cgf-gt-cgn") %>%
  filter(roi %in% hm1a.std.roi) %>%
  group_by(roi) %>%
  group_map(~ t.test(.$mean_topvoxels, 
                     alternative="greater",
                     mu = 0,
                     na.action("na.omit")))
names(hm1a.test.std) <- hm1a.std.roi

# Results Reported in Appendix E:
# "fear_lAmyg" ***
hm1a.test.std$fear_lAmyg
#M= 0.249536, t = 4.3344, df = 41, p-value = 4.623e-05

# "fear_lIFG"***
hm1a.test.std$fear_lIFG
#M=0.5513927, t = 5.2245, df = 41, p-value = 2.717e-06

# "fear_lThal" = NS
hm1a.test.std$fear_lThal 
#M=0.01899547, t = 0.45084, df = 41, p-value = 0.3272

# "fear_rAmyg" *
hm1a.test.std$fear_rAmyg
#M=0.1132042, t = 1.7425, df = 41, p-value = 0.04446

# "fear_rIFG" *
hm1a.test.std$fear_rIFG
#M=0.2187283, t = 1.981, df = 41, p-value = 0.02716

# "fear_rITG" = NS
hm1a.test.std$fear_rITG
#M=0.04033037, t = 0.44141, df = 41, p-value = 0.3306


### F > N (All Subjects), Robust One Sample T-Test - 1 ROI ####
hm1a.rob.roi <- c("fear_lmOccGyr")

hm1a.test.rob <- mc.fear %>%
  filter(contrast == "cgf-gt-cgn") %>%
  filter(roi %in% hm1a.rob.roi) %>%
  group_by(roi) %>%
  group_map(~ trdiff1sp(.$mean_topvoxels))

hm1a.test.rob  
#estimate (trimmed mean) =0.008268756
#test.stat 0.1045028
#p.value= 0.453





## Disgust Manipulation (H-M2a) ####

# Corresponds to pre-registered method (all ROIs normal) 

### D > N (All Subjects), One Sample T-Test - Reported for all ROI ####
hm2a.test <- mc.disg %>%
  filter(contrast == "cgd-gt-cgn") %>%
  group_by(roi) %>%
  group_map(~ t.test(.$mean_topvoxels, 
                     alternative="greater", 
                     mu = 0,
                     na.action("na.omit")))
names(hm2a.test) <- dspace

# Results: 
# "disgust_lAmyg" ***
hm2a.test$disgust_lAmyg
#M=0.5663547, t = 6.114, df = 41, p-value = 1.489e-07

# "disgust_lInsIFG" ***
hm2a.test$disgust_lInsIFG
#M=0.996387, t = 5.1524, df = 41, p-value = 3.43e-06

# "disgust_lPar" ***
hm2a.test$disgust_lPar
#M=1.502742, t = 8.2022, df = 41, p-value = 1.751e-10

# "disgust_rAmyg" ***
hm2a.test$disgust_rAmyg
#M=0.4162213, t = 4.3944, df = 41, p-value = 3.835e-05

# "disgust_rClaus" ***
hm2a.test$disgust_rClaus
#M=0.5057314, t = 5.4525, df = 41, p-value = 1.295e-06

# "disgust_rInsIFG" NS
hm2a.test$disgust_rInsIFG
#M=-0.2985867, t = -2.0419, df = 41, p-value = 0.9762




## Explicit vs. Implicit Responses (H-T1, H-P1, H-2) ####
### Add Covariates to Preference Intensity ####

table(s$partyID)

s.pref <- filter(s, match == 1) %>%
  select(participantID, Group, bathlaw,
         partyID, dem, rep, ind, 
         ideol, lib, cons, mod, 
         Self.ID, Safety, Jobs, Rights, 
         Id_culture, Moral, Health)
# Generate 7pt party scale to find leaners
s.pref$pid7 <- ifelse(s.pref$dem == 1, 1, 
                      ifelse(s.pref$dem == 2, 2,
                             ifelse(s.pref$rep == 1, 7,
                                    ifelse(s.pref$rep == 2, 6,
                                           ifelse(s.pref$ind == 1, 3,
                                                  ifelse(s.pref$ind == 2, 5,
                                                         ifelse(s.pref$ind == 3, 4, 99)))))))
#table(s.pref$pid7)
# Generate 7pt ideology scale to find leaners
s.pref$ideo7 <- ifelse(s.pref$lib == 1, 1, 
                       ifelse(s.pref$lib == 2, 2,
                              ifelse(s.pref$cons == 1, 7,
                                     ifelse(s.pref$cons == 2, 6,
                                            ifelse(s.pref$mod == 1, 3,
                                                   ifelse(s.pref$mod == 2, 5,
                                                          ifelse(s.pref$mod == 3, 4, 99)))))))
#table(s.pref$ideo7)
# Generate vars without (1) and with (2) leaners
s.pref$dem1 <- ifelse(s.pref$pid7 < 3, 1, 0)
s.pref$dem2 <- ifelse(s.pref$pid7 < 4, 1, 0)
s.pref$rep1 <- ifelse(s.pref$pid7 > 5 & s.pref$pid7 < 99, 1, 0)
s.pref$rep2 <- ifelse(s.pref$pid7 > 4 & s.pref$pid7 < 99, 1, 0)

s.pref$lib1 <- ifelse(s.pref$ideo7 < 3, 1, 0)
s.pref$lib2 <- ifelse(s.pref$ideo7 < 4, 1, 0)
s.pref$con1 <- ifelse(s.pref$ideo7 > 5, 1, 0)
s.pref$con2 <- ifelse(s.pref$ideo7 > 4, 1, 0)


#Generate variable for White self-ID
is.white <- c("White", "Caucasian", "White/Caucasian")
s.pref$white <- ifelse(s.pref$Self.ID %in% is.white, 1, 0)
#sum(s.pref$white)

# Generate binary of endorsement
s.pref$Endorse <- s.pref$Group - 1

### Add Binary Explicit Threat Ratings ####
# Make binary ratings, per pre-registered test (chi-sq)
s.pref$bin_safe <- ifelse(s.pref$Safety > 3, 1, 0)
s.pref$bin_jobs <- ifelse(s.pref$Jobs > 3, 1, 0)
s.pref$bin_rts <- ifelse(s.pref$Rights > 3, 1, 0)
s.pref$bin_id <- ifelse(s.pref$Id_culture > 3, 1, 0)
s.pref$bin_moral <- ifelse(s.pref$Moral > 3, 1, 0)
s.pref$bin_health <- ifelse(s.pref$Health > 3, 1, 0)

### Explicit Ratings Group Differences (H-T1 Tests) ####
# Pre-registration stipulates all contaminant threats (id, moral, health) are a family
# and that we will adjust critical thresholds for inference but 
# in practice this doesn't matter because none of the chi-sq tests is significant

# Contaminant threat 1 = id/culture
tab.id <- table(s.pref$Group, s.pref$bin_id)
set.seed(1234)
chisq.test(tab.id, simulate.p.value = T, B = 1000)
#X-squared = 4.4211, df = NA, p-value = 0.1239

# Contaminant threat 2 = morals
tab.moral <- table(s.pref$Group, s.pref$bin_moral)
set.seed(1234)
chisq.test(tab.moral, simulate.p.value = T, B = 1000)
#X-squared = 3.2308, df = NA, p-value = 0.2358

# Contaminant threat 3 = health
tab.health <- table(s.pref$Group, s.pref$bin_health)
set.seed(1234)
chisq.test(tab.health, simulate.p.value = T, B = 1000)
#X-squared = 1.0244, df = NA, p-value = 1

# For reference:
tab.safe <- table(s.pref$Group, s.pref$bin_safe)
set.seed(1234)
chisq.test(tab.safe, simulate.p.value = T, B = 1000)
#X-squared = 0, df = NA, p-value = 1

tab.jobs <- table(s.pref$Group, s.pref$bin_jobs)
set.seed(1234)
chisq.test(tab.jobs, simulate.p.value = T, B = 1000)
#X-squared = 1.0244, df = NA, p-value = 1

tab.rts <- table(s.pref$Group, s.pref$bin_rts)
set.seed(1234)
chisq.test(tab.rts, simulate.p.value = T, B = 1000)
#X-squared = 1.0244, df = NA, p-value = 1

# Exploratory: Social contamination index
sc_idx <- cbind(s.pref$bin_moral, s.pref$bin_id)
sc_idx <- rowSums(sc_idx)
s.pref$bin_scidx <- ifelse(sc_idx == 0, 0, 1)
tab.scidx <- table(s.pref$Group, s.pref$bin_scidx)
set.seed(1234)
chisq.test(tab.scidx, simulate.p.value = T, B = 1000)
# X-squared = 5.6757, df = NA, p-value = 0.05794

### Preference intensity (H-P1 Tests) ####
# Regions to model are those demonstrating between-group differences:
# lIFG (ER region), dmPFC (TomLoc region)
# Note that contrast values are scaled to help interpretation
# (mean of zero and sd of 1) 
# Most of the pre-registered analyses assumed distance measures would show group
# differences, but they did not

#### Add preference intensity to fmri data and subset to region and contrast ####
pref.m3.er <- left_join(x = m3.er, y = s.pref, by = c("participantID"))

d.lifg <- filter(pref.m3.er, 
                 roi == "er_lIFG",
                 contrast == "tgn-gt-cgn")

pref.m2 <- left_join(x = m2.tn, y = s.pref, by = c("participantID", "Group"))
d.dmpfc <- filter(pref.m2, 
                  ROI == "dmPFC",
                  extraction_contrast == "con_1_tgn-gt-cgn")

#### lIFG Models #####
# Baseline
ifg.base.s <- lm(bathlaw ~ scale(mean_topvoxels), d = d.lifg)
summary(ifg.base.s)

# Party ID, including leaners
ifg.pid.s <- update(ifg.base.s, . ~ . + dem2 + rep2)
summary(ifg.pid.s)

# Ideology, including leaners
ifg.ideo.s <- update(ifg.base.s, . ~ . + lib2 + con2)
summary(ifg.ideo.s)

# Race
ifg.race.s <- update(ifg.base.s, . ~ . + white)
summary(ifg.race.s)

ifg.full.s <- update(ifg.base.s, . ~ . +
                       dem2 + rep2 +
                       lib2 + con2 + white)
summary(ifg.full.s)

# Model Comparison for reference (small N so using Clarke Test) - Reported
# Note: requires clarkeTest package
clarkeTest::clarke_test(ifg.pid.s, ifg.ideo.s)
# Model 1 log-likelihood: -70
# Model 2 log-likelihood: -64
# Observations: 42
# Test statistic: 9 (21%)
# Model 2 is preferred (p = 0.00027)


#### dmPFC Models ####
# Baseline
dmpfc.base.s <- lm(bathlaw ~ scale(mean_extracted_betas), d = d.dmpfc)
summary(dmpfc.base.s)
# betas **

# Party ID, including leaners
dmpfc.pid.s <- update(dmpfc.base.s, . ~ . + dem2 + rep2)
summary(dmpfc.pid.s)

# Ideology, including leaners
dmpfc.ideo.s <- update(dmpfc.base.s, . ~ . + lib2 + con2)
summary(dmpfc.ideo.s)

# Race
dmpfc.race.s <- update(dmpfc.base.s, . ~ . + white)
summary(dmpfc.race.s)

# Full Model
dmpfc.full.s <- update(dmpfc.base.s, . ~ . + dem2 + rep2 +
                         lib2 + con2 + white)
summary(dmpfc.full.s)

### Information Equivalence (H-P2) - Threshold Criterion Not Met ####
# Our preregistration states that we will perform model comparison only if a 
# threat rating is a significant binary predictor of policy preference intensity.
# As this never occurs in the data (see below) we do not pursue this test further

# Id/Culture
mod0_id <- lm(bathlaw ~ bin_id, data = s.pref)
summary(mod0_id)

# Moral
mod0_moral <- lm(bathlaw ~ bin_moral, data = s.pref)
summary(mod0_moral)

# Health
mod0_health <- lm(bathlaw ~ bin_health, data = s.pref)
summary(mod0_health)



###################################

####################################
# 8. Figures (2a & 3a) ####
# Note: the replication data and files for figures 2b and 3b are 
# available in the Dataverse archive, but require special software (see README.txt)
## Fig 2a: M2 Violin Plot  ####
# Select the data for plotting
m2.plot <- m2.tn %>%
  select(participantID, Group, ROI, mean_extracted_betas) %>%
  data.frame()

# Format variables for plotting
m2.plot$Group <- factor(m2.plot$Group,
                        levels = c(1, 2),
                        labels = c("Oppose", "Endorse"))

# Create labels for significance test results
m2.plot.sig <- data.frame(ROI = unique(m2.plot$ROI),
                          Group = factor(c(rep(1, 7), rep(2,7)),
                                         levels = c(1, 2),
                                         labels = c("Oppose", "Endorse")),
                          Test = c("***", " ", "**", " ",
                                   " ", " ","",
                                   " "," "," ", " ",
                                   " "," "," "))
m2.plot.sig$Test <- as.character(m2.plot.sig$Test)

# Plot:
fig2a <- ggplot(data = m2.plot, 
                     aes(x = Group,
                         y = mean_extracted_betas,
                         group = Group,
                         fill = Group)) +
  geom_violin(alpha = 0.8) +
  geom_jitter(color = "black", 
              alpha=0.4,
              width = 0.05) +
  scale_fill_manual(values = c("white", "darkgrey")) +
  geom_hline(yintercept = 0) +
  labs(title="Activation in the Mentalizing Search Space",
       subtitle = "Points represent individual participants",
       y = "Mean Contrast (T>C)") +
  ylim(-2.8,4) +
  facet_wrap(.~ ROI,
             nrow = 2, ncol = 5) + 
  theme_bw() +
  theme(legend.position = "none",
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 16),
        strip.text = element_text(size = 14),
        title = element_text(size = 18))
fig2a
# Add stars 
fig2a_stars <- fig2a + 
  geom_text(data = m2.plot.sig, 
            mapping = aes(x = Group,
                          y = 3.8,
                          label = Test),
            size = 8) 
fig2a_stars


## Fig 3a: M3 Violin Plot ####

# Select the data for plotting
m3.plot <- m3.er %>%
  filter(contrast == "tgn-gt-cgn") %>%
  select(participantID, Group, roi, roi_labels, mean_topvoxels) %>%
  data.frame()
# Format variables for plotting
m3.plot$Group <- factor(m3.plot$Group,
                        levels = c(1, 2),
                        labels = c("Oppose", "Endorse"))
m3.plot$roi_labels <- factor(m3.plot$roi_labels,
                             labels = c("l. Cingulate","l. IFG",
                                        "l. MFG","l. MTG",
                                        "l. SFG","l. TG",
                                        "r. Insula/IFG","r. MFG",
                                        "r. Parietal","r. SFG"))
m3.plot$mean_topvoxels <- as.numeric(m3.plot$mean_topvoxels)

# Create labels for significance test results
m3.plot.sig <- data.frame(roi_labels = unique(m3.plot$roi_labels),
                          Group = factor(c(rep(1, 10), rep(2,10)),
                                         levels = c(1, 2),
                                         labels = c("Oppose", "Endorse")),
                          Test = c("*", "***", "*", "*","*",
                                   "*","*"," "," "," ",
                                   " "," "," "," "," ",
                                   " ","*"," "," "," "))
m3.plot.sig$Test <- as.character(m3.plot.sig$Test)

# Plot:
fig3a <- ggplot(data = m3.plot, 
                     aes(x = Group,
                         y = mean_topvoxels,
                         group = Group,
                         fill = Group)) +
  geom_violin(alpha = 0.8) +
  geom_jitter(color = "black", 
              alpha=0.4,
              width = 0.05) +
  scale_fill_manual(values = c("white", "darkgrey")) +
  geom_hline(yintercept = 0) +
  labs(title="Activation in the Self-Regulation Search Space",
       subtitle = "Points represent individual participants",
       y = "Mean Contrast (T>C)") +
  ylim(-2.8,4.2) +
  facet_wrap(.~ roi_labels,
             nrow = 2, ncol = 5) + 
  theme_bw() +
  theme(legend.position = "none",
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 16),
        strip.text = element_text(size = 14),
        title = element_text(size = 18))
fig3a

# Add stars option
fig3a_stars <- fig3a + 
  geom_text(data = m3.plot.sig, 
            mapping = aes(x = Group,
                          y = 3.8,
                          label = Test),
            size = 8) 
fig3a_stars 



#################################

#################################
# 9. Tables ####

## Appendix C: Study 1 Search Space Definition ####
# C1, C2: ToMLoc data and code available from authors on request

### Table C6: Fear Search Space Double Dissociations ####
froi.diss <- filter(fdn.diss, str_detect(roi, "fear"))

froi.df <- data.frame(ROI = c("l. Amygdala", "l. IFG", 
                              "l. Med. Occipital Gyrus", "l. Thalamus", 
                              "r. Amygdala", "r. IFG",
                              "r. ITG"),
                      FN_m = round(froi.diss$FN_trm_tr1s, 2),
                      FN_t = round(froi.diss$FN_stat_tr1s, 2),
                      FN_p = round(froi.diss$FN_pvalue_tr1s,2),
                      FD_m = round(froi.diss$FD_trm_tr1s, 2),
                      FD_t = round(froi.diss$FD_stat_tr1s, 2),
                      FD_p = round(froi.diss$FD_pvalue_tr1s,2))
froi.df

### Table C7: Disgust Search Space Double Dissociations ####
# Uses fdn.diss from above (also saved out as csv)
droi.diss <- filter(fdn.diss, str_detect(roi, "disgust"))

droi.df <- data.frame(ROI = c("l. Amygdala", "l. Insula/IFG", 
                              "l. Parietal", "r. Amygdala", 
                              "r. Claustrum", "r. Insula/IFG"),
                      DN_m = round(droi.diss$DN_trm_tr1s, 2),
                      DN_t = round(droi.diss$DN_stat_tr1s, 2),
                      DN_p = round(droi.diss$DN_pvalue_tr1s, 2),
                      DF_m = round(droi.diss$FD_trm_tr1s, 2),
                      DF_t = round(droi.diss$FD_stat_tr1s, 2),
                      DF_p = round(droi.diss$FD_pvalue_tr1s, 2))
droi.df

## Appendix D: Study 1 Additional Tables ####
### M1: Disgust and Fear T Condition Similarities ####

#### Table D1: Endorsers, Disgust ROIs ####
tab_d1.df <- data.frame(ROI = c("l. Amygdala", "l. Insula/IFG", 
                                "l. Parietal", "r. Amygdala", 
                                "r. Claustrum", "r. Insula/IFG"),
                        TDTN_diff = round(h1a.m1.disgust$TDTN_diff, 2),
                        TDTN_p = round(h1a.m1.disgust$TNTD_pvalue, 2),
                        TDTF_diff = round(h1a.m1.disgust$TDTF_diff, 2),
                        TDTF_p = round(h1a.m1.disgust$TDTF_pvalue, 2),
                        H1a_disgust = h1a.m1.disgust$D_sim)
tab_d1.df

#### Table D2: Endorsers, Fear ROIs ####
tab_d2.df <- data.frame(ROI = c("l. med. Occipital Gyrus","r. Amygdala", "r. ITG"),
                        TFTN_diff = round(h1a.m1.fear$TFTN_diff, 2),
                        TFTN_p = round(h1a.m1.fear$TNTF_pvalue, 2),
                        TFTD_diff = round(h1a.m1.fear$TFTD_diff, 2),
                        TFTD_p = round(h1a.m1.fear$TDTF_pvalue, 2),
                        H1a_fear = h1a.m1.fear$F_sim)
tab_d2.df

#### Table D3: Opponents, Disgust ROIs ####
# Uses 0.2 trim to match g2)
tab_d3.df <- data.frame(ROI = c("l. Amygdala", "l. Insula/IFG", 
                                "l. Parietal", "r. Amygdala", 
                                "r. Claustrum", "r. Insula/IFG"),
                        TNTD_diff = round(h1b.m1.disgust$TNTD_diff, 2),
                        TNTD_p = round(h1b.m1.disgust$TNTD_pvalue, 2),
                        TNTF_diff = round(h1b.m1.disgust$TNTF_diff, 2),
                        TNTF_p = round(h1b.m1.disgust$TNTF_pvalue, 2),
                        H1b_test = h1b.m1.disgust$N_sim)
tab_d3.df


#### Table D4: Opponents, Fear ROIs ####
tab_d4.df <- data.frame(ROI = c("l. Amygdala","l. IFG",
                                "l. med. Occipital Gyrus",
                                "l. Thalamus",
                                "r. Amygdala", "r. IFG","r. ITG"),
                        TNTF_diff = round(h1b.m1.fear$TNTF_diff, 3),
                        TNTF_p = round(h1b.m1.fear$TNTF_pvalue, 3),
                        TNTD_diff = round(h1b.m1.fear$TNTD_diff, 3),
                        TNTD_p = round(h1b.m1.fear$TNTD_pvalue, 3),
                        H1b_test = h1b.m1.fear$N_sim)
tab_d4.df
#write.csv(tab_d4.df, "Opp_Fear_TSims_final.csv", row.names=F)


#### Table D5: Group Differences, Disgust ROIs ####
tab_d5.df <- data.frame(ROI = c("l. Amygdala", "l. Insula/IFG", 
                                 "l. Parietal", "r. Amygdala", 
                                 "r. Claustrum", "r. Insula/IFG"),
                         G1_tm = round(h1c.m1.df$g1_tm[1:6], 2),
                         G2_tm = round(h1c.m1.df$g2_tm[1:6], 2),
                         y_stat = round(h1c.m1.df$y_stat[1:6], 2),
                         y_pval = round(h1c.m1.df$y_pval[1:6], 2))
tab_d5.df


#### Table D6: Group Differences, Fear ROIs ####
tab_d6.df <- data.frame(ROI = c("l. Amygdala", "l. IFG", 
                              "l. Med. Occipital", "l. Thalamus", 
                              "r. Amygdala", "r. IFG",
                              "r. ITG"),
                         G1_tm = round(h1c.m1.df$g1_tm[7:13], 2),
                         G2_tm = round(h1c.m1.df$g2_tm[7:13], 2),
                         y_stat = round(h1c.m1.df$y_stat[7:13], 2),
                         y_pval = round(h1c.m1.df$y_pval[7:13], 2))
tab_d6.df

### M2: Mentalizing ####
#### Table D7: Opponents ####
tab_d7.df <- data.frame(ROI = h2a.test.df$ROI,
                        TN_Contrast = round(h2a.test.df$estimate, 3),
                        t_statistic = round(h2a.test.df$statistic, 3),
                        p_value = round(h2a.test.df$p.value, 3))
tab_d7.df

#### Table D8: Endorsers ####  
tab_d8.df <- data.frame(ROI = h2a.check.df$ROI,
                        TN_Contrast = round(h2a.check.df$estimate, 3),
                        t_statistic = round(h2a.check.df$statistic, 3),
                        p_value = round(h2a.check.df$p.value, 3))
tab_d8.df

#### Table D9: Group Differences ####
tab_d9.df <- data.frame(ROI = h2b.tests.df$ROI,
                        G1_Contrast = round(h2b.tests.df$mean.of.x, 2),
                        G2_Contrast = round(h2b.tests.df$mean.of.y, 2),
                        t_statistic = round(h2b.tests.df$statistic, 2),
                        p_value = round(h2b.tests.df$p.value, 2))
tab_d9.df

### M3: Self-Reg ####
#### Table D10: Opponents ####

tab_d10.df <- data.frame(ROI = h3a.test.df$ROI,
                        TN_Contrast = round(h3a.test.df$estimate, 3),
                        t_statistic = round(h3a.test.df$t.stat, 3),
                        p_value = round(h3a.test.df$p.value, 3))
tab_d10.df



#### Table D11: Endorsers ####
tab_d11.df <- data.frame(ROI = h3a.check.df$ROI,
                         TN_Contrast = round(h3a.check.df$estimate, 3),
                         t_statistic = round(h3a.check.df$t.stat, 3),
                         p_value = round(h3a.check.df$p.value, 3))
tab_d11.df



#### Table D12: Group Differences ####
tab_d12.df <- data.frame(ROI = h3b.test.df$ROI,
                         G1_TN_Contrast = round(h3b.test.df$g1.mean, 3),
                         G2_TN_Contrast = round(h3b.test.df$g2.mean, 3),
                         t_statistic = round(h3b.test.df$t.stat, 3),
                         p_value = round(h3b.test.df$p.value, 3))
tab_d12.df



### Additional Hypotheses ####
#### Table D13: left IFG Regression Models ####

# N= 42 
mod1.univ <- ifg.base.s
mod2.univ <- update(ifg.base.s, . ~ . + dem2 + rep2) #incl. leaners
mod3.univ <- update(ifg.base.s, . ~ . + lib2 + con2) #incl. leaners
mod4.univ <- update(ifg.base.s, . ~ . + white)
mod5.univ <- update(ifg.base.s, . ~ . + dem2 + rep2 + 
                      lib2 + con2 + white)
summary(mod5.univ)


vars.univ <- c("Mean Contrast Value (T>C)", "Democrat", "Republican",
               "Liberal", "Conservative", "White", "Constant")

stargazer::stargazer(
  mod1.univ, mod2.univ, mod3.univ, mod4.univ, mod5.univ,
  type = "text",
  covariate.labels = vars.univ,
  dep.var.caption = c("Bathroom Bill Support (1-5)"),
  dep.var.labels = c("Strongly Oppose=1, Strongly Endorse=5"),
  column.labels = c(rep("",5)),
  star.cutoffs = c(0.05, 0.01, 0.001),
  digits = 2,
  df = F,
  omit.stat = c( "ser"),
  column.sep.width = "-0pt" )


#### Table D14: DMPFC Regression Models ####

# N= 41
mod1.univ <- dmpfc.base.s
mod2.univ <- update(dmpfc.base.s, . ~ . + dem2 + rep2) #incl. leaners
mod3.univ <- update(dmpfc.base.s, . ~ . + lib2 + con2) #incl. leaners
mod4.univ <- update(dmpfc.base.s, . ~ . + white)
mod5.univ <- update(dmpfc.base.s, . ~ . + dem2 + rep2 + 
                      lib2 + con2 + white)
summary(mod5.univ)


vars.univ <- c("Mean Contrast Value (T>C)", "Democrat", "Republican",
               "Liberal", "Conservative", "White", "Constant")

stargazer::stargazer(
  mod1.univ, mod2.univ, mod3.univ, mod4.univ, mod5.univ,
  type = "text",
  covariate.labels = vars.univ,
  dep.var.caption = c("Bathroom Bill Support (1-5)"),
  dep.var.labels = c("Strongly Oppose=1, Strongly Endorse=5"),
  column.labels = c(rep("",5)),
  star.cutoffs = c(0.05, 0.01, 0.001),
  digits = 2,
  df = F,
  omit.stat = c( "ser"),
  column.sep.width = "-0pt" )

## Appendix F: Study 2 Additional Materials ####
# Note: Study 2 raw data is only available from the original authors (see Study 2 RepCode for details)
# This code reproduces statistics reported in Columns 3 and 4 of Table F1 

### Table F1: Study Demographics (Study 1 only) ####
s.match <- filter(s, match == 1)

# Gender (same for both groups)
round(table(s.match$Group, s.match$Sex)[1,]/21,2)
round(table(s.match$Group, s.match$Sex)[2,]/21,2)

# Age
mean(filter(s.match, Group == 1)$age)
sd(filter(s.match, Group == 1)$age)
mean(filter(s.match, Group == 2)$age)
sd(filter(s.match, Group == 2)$age)

# White (Caucasian + White + White/Caucasian)
w.selfid <- c("Caucasian", "White", "White/Caucasian")
w.subs <- filter(s.match, Self.ID %in% w.selfid)
table(w.subs$Group)/21

# 4-category Party ID (1 = Democrat, 2 = Republican, 3 = Independent, 4 = DK/Other)
table(s.match$Group, s.match$partyID)[1,]/21
table(s.match$Group, s.match$partyID)[2,]/21
# Note that although regression analyses include leaners,
# this table reports without leaners to be comparable to 
# NASIS data, which does not ask 2-part PID

# 4-category Ideology (1 = Liberal, 2 = Conservative, 3 = Moderate, 4 = DK/Other)
table(s.match$Group, s.match$ideol)[1,]/21
table(s.match$Group, s.match$ideol)[2,]/21
# Same caveat as party ID

